A mechanical model of ocular bulb vibrations and implications for acoustic tonometry

In this study, we propose a comprehensive mechanical model of ocular bulb vibrations and discuss its implications for acoustic tonometry. The model describes the eye wall as a spherical, pre-stressed elastic shell containing a viscoelastic material and accounts for the interaction between the elastic corneoscleral shell and the viscoelastic vitreous humor. We investigate the natural frequencies of the system and the corresponding vibration modes, expanding the solution in terms of scalar and vector spherical harmonics. From a quantitative point of view, our findings reveal that the eyebulb vibration frequencies significantly depend on IOP. This dependency has two origins: “geometric” stiffening, due to an increase of the pre-stress, and “material” stiffening, due to the nonlinearity of the stress-strain curve of the sclera. The model shows that the second effect is by far dominant. We also find that the oscillation frequencies depend on ocular rigidity, but this dependency is important only at relatively large values of IOP. Thus close to physiological conditions, IOP is the main determinant of ocular vibration frequencies. The vitreous rheological properties are found to mostly influence vibration damping. This study contributes to the understanding of the mechanical behavior of the eye under dynamic conditions and thus has implications for non-contact intraocular pressure measurement techniques, such as acoustic tonometry. The model can also be relevant for other ocular pathological conditions, such as traumatic retinal detachment, which are believed to be influenced by the dynamic behavior of the eye.


Motion equations
We model the eye as a composite structure, consisting of a viscoelastic material, the vitreous humor, which occupies a spherical domain B R of radius R centered at the origin, enclosed within a thin elastic shell, representing the corneo-scleral complex.We consider an equilibrium state whereby the vitreous humor is at a pressure p, sustained by a state of stress in the corneo-scleral shell, and we study small-amplitude vibrations about this state.
The equations that govern the evolution of the displacement u in the vitreous humor are where ρ v is the mass density of the vitreous humor, a superposed dot denotes partial differentiation with respect to time, and S is the increment of the Piola (nominal) stress with respect to the reference state.For the incremental stress we adopt the following constitutive equations: where q is the pressure increment, and D( u) = sym dev ∇ u is the symmetric deviatoric part of the velocity gradient.
In the first of (2), the second term on the right-hand side accounts for the change of the Piola stress due to the initial pressure.This term, whose derivation may be found in [2], does not make any contribution to the bulk equations, since div∇u ⊤ = ∇div u = 0; however, it makes a contribution to the boundary conditions, as shown below.The relaxation function G(t) describes the mechanical behavior of the material in the linear setting.For an elastic material, the relaxation function is proportional to the Heaviside function.For a viscous material, G(t) is the Dirac delta centered at the origin.
We assume that there is no slip between the vitreous humor and the sclera.Accordingly, we set the displacement field on the sclera to be equal to the vector field u on the sphere S R = ∂B R .For convenience, we decompose u into its normal and tangential components, respectively, w and v. Accordingly, on denoting by n the unit normal to S R , we write The equations that govern small displacements of the shell from its reference configuration have been derived in [5], and are given by Here ρ s is the density of the sclera, h is its thickness, P = I − n ⊗ n is the projector on the plane orthogonal to n, and b = −Sn, is the increment of the nominal traction (i.e., the traction per unit reference area) acting on the shell.By (4) and ( 5) we have The tensor fields N and M represent the increment of the nominal membrane-force and bending-moment tensors, and they are given by the constitutive equations where µ and  λ are suitable material moduli, which depend on the nature of the material comprising the shell, and respectively, the stretching and bending strains.The symbol ∇ s denotes the surface gradient on S R .The surface divergence of the tangential vector field v is defined as div s v = tr(∇ s v) = P : ∇ s v, where the symbol : denotes contraction between tensors, while the surface divergence of N is defined by requiring that a • div s N = div s (N ⊤ a) for every constant vector field a (note that since Nn = 0, N ⊤ a is a tangential vector field).The same definition holds for div s M.
It is worth remarking that the first term on the right-hand side of (7) and the first term on the right-hand side of (2) share the same nature: they take into account increments of the nominal stress that originate from the reference state being pre-stressed.
We next elaborate the last term in each of the motion equations (6).First we compute since ∇n = 1 R P. Thus, where ∇ s w is the superficial gradient of w.Accordingly, we have from ( 2) Then, the motion equations (6) can be rewritten as Harmonic solutions Frequency domain.We look for solutions of (1), (2), of the form For any such solution, we have S(x, t) = Re(e ζt S(x, ζ)), where S = p∇u ⊤ + Σ, with and Thus, system (1) takes the form: By substituting (7), (8), and (11) into (12), keeping (13) in mind: where w = u • n and v = Pu, are, respectively, the normal and the tangential components of u.
Expansion in spherical harmonics.
To solve (16), we mimick the procedure of [3].For r the position vector with respect to the center of B R , we let r = |r| and  r = r/r, and we introduce the expansion of the pressure increment using spherical harmonics {Y ℓm : ℓ ≥ 0, −ℓ ≤ m ≤ ℓ}.Spherical harmonics are a complete orthonormal set on the unit sphere.Using spherical polar coordinates (θ, φ), with 0 ≤ θ ≤ π the polar angle (θ = 0 on the North Pole) and 0 ≤ φ < 2π the azimuthal angle, spherical harmonics are defined by where N ℓm is a normalization constant and P ℓm : [−1, +1] → R is an associated Legendre polynomial.For m ≥ 0, the associated Legendre polynomials are defined by , with P ℓ (x) the ℓ-th Legendre polynomial, defined by P 0 (x) = 1, P 1 (x) = x, and by the recursion formula (ℓ + 1) (l+m)! P ℓm (x) [4].The analogous of spherical harmonics for vector fields are the vector spherical harmonics, defined on the unit sphere through [4, p. 1898]: where ∇ s is the surface gradient on the unit sphere, and We remark that the first spherical harmonic Y 00 ( r) is constant and, as a consequence, Vector spherical harmonics are a complete orthonormal set on the unit sphere.Accordingly, the displacement u(r, ζ) admits the expansion: with where, in view of (21), we have set Differential identities for spherical harmonics.With a view towards solving (16), we extend spherical harmonics and vector spherical harmonics to the entire space by letting Y ℓm (r) = Y ℓm ( r), p ℓm (r) = p ℓm ( r), etc.The resulting extensions satisfy [1] ∇Y Moreover, and

Solution in the vitreous
The identities (25)-( 27) make the choice of spherical harmonics particularly useful to solve (16)-( 17).In fact, by ( 25)-( 27), we have where Furthermore, by making use of (26), we find Moreover, on setting We remark that, although Y 00 is constant, the vector spherical harmonic p 00 (r) = Y 00  r is not.Now, the key observation is that the action on the spherical harmonics of the operators ∆, div, and ∇, which appear in the motion equations (16), does not mix coefficients with different values of ℓ and m.Keeping this observation in mind, we substitute the expansions (19) and ( 22) into the governing equations (16) and we use (28) and (31), as well as the aforementioned orthogonality property to obtain a series of uncoupled systems, each for a particular value of ℓ ≥ 0 and −ℓ ≤ m ≤ ℓ.
The resulting systems have to be treated differently, according to whether ℓ = 0 or ℓ ≥ 1.As we shall see below, the corresponding solutions are trivial in the first case.
Coefficients with ℓ=0.For ℓ = 0 we obtain the following system for P 00 (r, ζ) and Q 00 (r, ζ): All solutions of the second equation are proportional to 1/r 2 .Thus, the requirement that solutions are bounded demands that P 00 (r, ζ) = 0.Then, the first of (34) imposes that Q 00 (r, ζ) be a constant.Thus, (1) 00 an integration constant.This solution describes the equilibrium configuration with an arbitrary pressure.
The case ℓ ≥ 1.For each ℓ ≥ 1 and −ℓ ≤ m ≤ ℓ, we multiply both sides of the first of (16) by, respectively, p ℓm , b ℓm and b ℓm , and we multiply both sides of the second of ( 16) by Y ℓm .Integrating on the unit sphere and making use of the aforementioned orthogonality properties, we obtain a system of four differential equations: where we have set The first three equations of (34) are uncoupled from the last one; their solution can be obtained by elimination of the unknowns B ℓm (r, ζ) (using the third equation) and Q ℓm (r, ζ) (using the second equation).This leads to the fourth-order differential equation for P ℓm (r, ζ): The solution of (36) depends, in general, on four integration constants.These constants determine B ℓm (r, ζ) and Q ℓm (r, ζ), respectively through the second and the third of (34).Among these solutions, we retain as admissible those which are bounded.Likewise, the fourth equation of (34) determines C ℓm (r, ζ) to within two integration constants.The result is: where j ℓ (x) =  π/(2x)J ℓ+1/2 (x) denotes the ℓ-th spherical Bessel function, (1) ℓm and C (3) ℓm are three constants.Concerning the coefficients with ℓ = 1, the smoothness of the displacement at the origin imposes (see [3] and references cited therein) that B 1m (0) = √ 2P 1m (0) and P ′ 1m (0) = B ′ 1m (0) = 0.The first two spherical Bessel functions, for small values of the argument z, can be approximated as Thus, near r = 0, (40) 1m and C (2) 1m .On the other hand, the requirement that B 1m (0) = √ 2P 1m (0) yields the condition As a result, we have and using the properties of the spherical harmonics Y 1m one can check that the resulting displacement is uniform (i.e., a translation).As one can guess, the only oscillatory motion compatible with a uniform translation is the trivial one, and this shall be confirmed in the analysis carried out in the next section.

Matching conditions.
The integration constants C (i) ℓm are determined by enforcing the matching conditions (17).To convert (17) into conditions at r = R on the functions Q ℓm (r, ζ), P ℓm (r, ζ), B ℓm (r, ζ), and C ℓm (r, ζ), we proceed in two steps.In the first step, we use ( 19) and ( 22) to compute the left-hand sides of the matching conditions (16).The calculation is quite straightforward, and leads to: and where we have set b 0m = 0 and c 0m = 0.
In the second step, we substitute the expansions (22) into the right-hand sides of the matching conditions (16).This step is a bit more elaborate than the previous one, because there are numerous terms to compute.To make the calculation more organized, we write the tangential and normal components of the displacement on S R as where and Then, we compute the action on v and w of the linear operators appearing on the right-hand sides of ( 16) by superposition on each term of the summations in (45): To deduce (48), it suffices to notice that the identities (25)-( 27) hold also when the differential operators ∇, div and ∆ are replaced by the surface differential operators ∇ s , div s and ∆ s acting on the restrictions of their argument on the sphere S R .The calculations are routine, but for completeness, we prefer to work them out.Since Y ℓm does not depend on the radial distance r from the origin, on S R we have ∇Y ℓm • n = 0, and therefore, Moreover, on denoting by ∂ r the derivative in the radial direction, we have, on S R , Since ∇ s Y ℓm is orthogonal to n, we have ∂ r ∇ s Y ℓm •n = 0, thus (recall that : is the contraction between tensors): In a similar fashion, since ∂ r p ℓm = 0, divp ℓm = tr(∇p ℓm ) = tr(∇ s p ℓm ) =: div s p ℓm .
It is also easy to check that divb ℓm = div s b ℓm , divc ℓm = div s c ℓm = 0. To show that ∆p ℓm = ∆ s p ℓm , we take a constant vector a and we apply the definition of divergence of a tensor field to ∇p ℓm : where we use the fact that b ℓm does not depend on r.Moreover, (50) Analogous calculations hold for b ℓm and c ℓm .With the above observation, and resorting to the identities (25)-( 27) we can compute the action on spherical harmonics of the operators appearing on the right-hand sides of the boundary conditions (17): Substituting (45) into the right-hand sides of (17) we obtain (53) Here, again, we have adopted the convention that b 0m = 0 and c 0m = 0.
Equating the right-hand sides of ( 43)-( 44) with the right-hand sides of ( 52)-( 53) and using the orthogonality properties of spherical and vector spherical harmonics we obtain, for ℓ = 0, the equation and, for ℓ ≥ 1, the following system: where the coefficients are given by Note that the left-hand sides of (55) have physical dimensions of pressure/time.

Characteristic equation
We now substitue into (54) and (55) the representation formulas (33) and (37).Again, we must treat separately the cases ℓ = 0, and ℓ ≥ 1.The case ℓ = 0 is however trivial.In fact, in view of (33), (54) yields C (1) 00 = 0. Thus we conclude that P 00 (r, ζ) = 0 and Q 00 (r, ζ) = 0.For ℓ ≥ 1 we substitute the representation (37) into (55), and we use the recurrence properties of the spherical Bessel functions to the homogeneous system where and the homogeneous equation where For ℓ = 1, the system (57) is supplemented by the regularity condition (41), which we rewrite here for the reader's convenience: Clearly, (57) and (61) constitute an overdetermined system, which admits only the trivial solutions, consistent with the expectation that during an oscillatory process the motion of the system cannot be a uniform displacement.Thus, the equation for ℓ ≥ 2, and the equation for ℓ ≥ 1 select the values of ζ for which oscillatory motions are possible.Summing up, we obtain two family of modes.The modes of the first family, for ℓ ≥ 2, are a linear combination of the vector spherical harmonics p ℓm and b ℓm , which possess both a radial and a tangential component.Modes of the second family, for ℓ ≥ 1, are proportional to the vector spherical harmonics c ℓm , which are purely tangential.We observe that the mode of the second family corresponding ℓ = 1 does not entail a deformation of the shell.
The case of a pressurized spherical shell.The numerical solution of (63) and (62) involves a root-finding algorithm that requires an initial guess.To obtain this guess, we use the vibration frequencies from the uncoupled problem where a pressurized spherical shell vibrates freely.To obtain the equations that give these frequencies, we set to zero the right-hand sides of (17), which account for the effect of the material filling the shell.In terms of spherical harmonics, this is equivalent to setting to null the right-hand sides of ( 52) and (53).Then, on projecting these equations on each spherical harmonic, we obtain the following homogeneous system: The first is a quadratic equation with respect to the variable ζ 2 .Writing the explicit form of the first equation would require too much space, and is not particularly informative.However, it is useful to note that for ℓ = 1, the first equation reduces to It is immediate that they vanish for ℓ = 1, since in that case s 2 ℓ = 2.The corresponding solution describes rigid rotations with constant angular velocity.